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ABSTRACT 

The  computer  simulation  methods  have  been  applied  to  study  the  structure  of  aqueous  solutions  of  simple  ionic  salts  in 
the  region  of  very  high  concentrations.  The  calculations  of  ionic  structures  in  solutions  were  performed  for  NaOH,  NaCl, 
LiCl  and  MgCl2  solutions.  The  concentrations  ranged  from  0.2  M to  saturated  solutions,  in  some  cases  as  much  as  19  M. 

A particularly  careful  analysis  was  devoted  to  the  topology  of  the  ionic  structures  in  solution.  Up  to  now,  most  of  the 
research  on  ionic  solutions  was  devoted  to  studies  of  the  ionic  hydration  shells.  However,  beyond  the  Debye  - Huckel  range 
of  very  low  concentrations,  very  little  is  known  about  the  interionic  spatial  correlations.  Certain  theories  predict  the 
existence  of  ionic  quasi-lattices  in  the  region  of  high  concentrations. 

In  the  present  work  we  used  the  Molecular  Dynamics  method  combined  with  such  statistical  tools  as  the  radial 
distribution  functions,  Voronoi  tessellations,  the  running  and  O'Keeffe  coordination  numbers,  etc.,  in  order  to  analyze  the 
ionic  structures.  The  radial  distribution  functions  of  three  types:  the  cation  - anion,  cation  - cation  and  the  anion  - anion 
type  were  calculated  for  each  solution.  The  functions  are  typical  for  the  quasi-crystalline  order  within  the  first  2-3  ionic 
coordination  layers  around  a selected  ion.  The  order  is  particularly  pronounced  for  the  anion  - cation  RDF’s. 

The  distributions  of  the  sphericity  factor  of  the  Voronoi  polyhedra  were  calculated  for  the  ionic  substructures  in  the 
configurations  produced  by  the  Molecular  Dynamics  simulation.  The  increase  of  the  ion  concentration  causes  evolution  of 
these  distributions  towards  increased  signatures  of  predominant  geometries  of  the  Voronoi  polyhedra.  This,  together  with 
the  results  for  RDF’s,  provides  a strong  conjecture  for  existence  of  the  ordered  structures  of  ions  in  concentrated  solutions. 


Keywords:  ionic  solutions,  lattice  theory,  computer  simulation. 


1.  INTRODUCTION 

The  structure  of  ionic  solutions  - both  the  local  structure  of  the  solvation  layers  around  the  ions  and  the  topology  of  the 
ion  distributions  in  solution  - is  of  primary  importance  for  our  understanding  of  nucleation  phenomena  and  growth  of 
crystals  from  the  solutions.  Unfortunately  the  aqueous  solutions  of  electrolytes  represent  a complicated  case  and  can  not  be 
described  by  any  analytical  theory  in  the  full  range  of  concentrations.  For  low  concentrations  of  ions  there  is  a number  of 
theoretical  models  that  work  reasonably  well:  among  others  the  Debye-Htickel  theory,  the  hyper-netted  chain  theory,  the 
mean  spherical  theory  (see  u 2 for  references).  All  these  theories  converge  to  the  same  result  - the  limiting  laws  of  the 
Debye-Htickel  theory  - when  the  ion  concentration  decreases  to  zero.  On  the  other  hand,  at  the  range  of  high  concentrations 
the  validity  of  the  models  deteriorates  very  quickly  and  another  theory  should  be  introduced.  Such  a theory,  based  on  the 
lattice-like  topology  of  the  ionic  structures,  has  been  mentioned  and  discussed  many  times  in  the  literature.  The  consequent 
formulation  of  the  theory  has  been  provided  by  Frank  and  Thompson  3.  The  most  complete  version  of  the  lattice  theory  of 
concentrated  electrolytes  was  given  in  a series  of  papers  by  Ruff  and  coworkers 4’ 5’ 6.  The  model  proposed  by  Ruff  is  can  be 
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described  as  follows.  (1)  The  ions  are  distributed  in  a lattice  - like  arrangement.  The  positions  of  the  ions  are  allowed  to 
deviate  randomly  from  the  lattice  sites.  (2)  The  ions  are  immersed  in  a structureless  continuous  dielectric  and  their 
coulombic  interactions  are  scaled  down  with  the  average  dielectric  constant  which  may  be  taken  as  a function  of  interionic 
separation  distance.  The  medium  around  the  ions  has  a dielectric  gradient  which  acts  as  a repulsive  force  between  the  ions  at 
short  separation  distances. 

Comparisons  of  calculated  and  observed  data  for  the  dielectric  constant,  activity  coefficients  and  partial  molar 
enthalpy  of  several  1-1  and  1-2  electrolytes  show  fairly  good  agreement  in  entire  concentration  range  4,s.  It  has  been  shown  6 
that  in  the  limit  of  dilute  systems,  the  deviations  of  the  ion  distribution  from  the  lattice  positions  grow  and  the  lattice  theory 
converges  to  the  Debye-Hiickel  theory. 

In  die  modem  investigations  of  liquids  and  liquid  solutions  the  computer  experiments  play  the  important  role 
occupying  in  a way  an  intermediate  place  between  theory  and  experiment.  In  the  case  of  very  complicated  problems,  that 
are  impossible  to  treat  in  an  analytical  way,  the  simulations  correspond  to  the  theory.  On  the  other  hand,  for  simpler  and 
more  tractable  problems,  the  simulations  can  be  used  as  a check  against  existing  analytical  theories.  Hence,  we  propose  to 
check  the  lattice  theory  of  concentrated  strong  electrolytes  by  comparing  it  with  the  results  of  the  computer  experiment  for 
these  media. 

The  computer  simulation  method  has  been  frequently  used,  with  remarkable  success,  to  predict  the  structure  and 
physico-chemical  properties  of  aqueous  electrolyte  solutions.  The  first  simulations  were  performed  by  Heinzinger  as  early 
as  1974.  We  refer  to  the  review  papers 7 8 as  the  source  of  references.  However,  with  very  few  exceptions  like  e.g.  9,  most  of 
the  calculations  concerned  the  solvation  structures  of  the  ions,  that  is  the  geometrical  arrangements  of  water  molecules  in 
the  vicinity  of  the  cations  and  anions.  The  ion  - ion  correlations  were  either  neglected  or  shown  by  the  radial  correlation 
functions.  The  investigation  of  the  geometry  of  the  ion  arrangements  is  very  difficult.  In  order  to  obtain  reliable  results  one 
has  to  include  in  the  simulation  box  at  least  a hundred  of  ions  and  proportional  number  of  water  molecules  that  results  from 
relative  concentrations.  This  increases  the  quantity  of  atomic  objects  in  the  simulation  box  to  almost  intractable  number  and 
such  calculations  have  become  possible  only  recently  due  to  the  increase  of  the  speed  of  available  computers. 


2.  COMPUTER  SIMULATION  METHOD 


2.1.  Simulation  details. 

The  simulations  have  been  performed  by  the  Molecular  Dynamics  method  10  using  the  flexible  model  of  water.  The 
box  contained  400-2000  water  molecules  and  the  number  of  ions  that  results  from  the  concentration  of  the  solution  and  the 
size  of  the  box  (16-160  ions).  The  simulation  box  was  assumed  in  the  cubic  form  with  the  size  ranging  between  2.5  and  4.1 
nm,  depending  on  the  ion  concentration.  The  calculations  were  performed  for  NaCl  and  MgCl2  solutions  with 
concentrations  ranging  from  0.2  to  about  5.0M  as  well  as  for  LiCl  and  NaOH  solutions  with  concentrations  from  0.5  to 
19.0M. 

The  ion-ion  interaction  was  assumed  as  a sum  of  the  Coulomb  potential  and  the  Lennard-Jones  interactions.  The 
interaction  parameters  were  taken  from  the  Universal  Force  Field  version  1.02  and  the  shifted  force  method  for 
non-coulombic  interactions  has  been  applied.  The  periodic  boundary  conditions  have  been  employed  with  the  Ewald 
summation  for  calculations  of  the  long-range  electrostatic  forces.  A typical  simulation  run  included  three  stages:  (a)  the 
initial,  equilibration  phase  of  l+2xl04  time  steps  with  the  temperature  scaling,  (b)  the  control  phase  of  lxlO4  time  steps  and 
(c)  the  proper  simulation  of  4+8.5xl04  time  steps.  The  elapsed  time  of  the  proper  simulation  reached  40+85  ps  and  the 
temperature  oscillated  near  300  K.  The  analysis  of  the  ionic  structures  was  performed  using  the  radial  distribution  functions 
(RDF)  and  stochastic  Voronoi  polyhedra. 

2.2.  Analysis  of  radial  distribution  functions  and  Voronoi  polyhedra. 

The  local  structure  of  a disordered  medium  can  be  described  by  a set  of  distribution  functions  for  atomic  positions.  The 
pair  distribution  function  gifarj)  (usualy  calculated  as  a function  g(r)  of  radial  distance  only  and  called  the  radial 
distribution  function,  RDF)  is  most  often  used.  The  value  of  g(r)  at  a given  r represents  the  density  of  probability  of  finding 
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a pair  of  atoms  separated  by  the  distance  r relative  to  the  probability  expected  for  an  uniform  distribution  of  atoms  with  the 
same  density.  The  formula  for  calculating  g(r)  is  given  by: 


g(r)  = (V/N)2  < SZy  m 5(r  - nj)  > 


(1) 


where  N is  the  number  of  particles  in  volume  V,  8 denotes  the  Dirac  delta,  the  brackets  < > denote  the  averaging  over  the 
ensemble  of  the  particles,  the  sums  over  / and  j run  over  the  different  pairs  of  the  particles  in  the  simulation  box. 

For  a multicomponent  system  the  radial  distribution  function  can  be  defined  for  each  type  of  the  atom  pair.  For  an  AB 
mixture  the  following  functions  : gM,  gM,  gBs  can  be  defined.  In  the  present  paper  we  calculate  the  RDF’s  describing 
spatial  correlations  between  ions:  the  cation  - anion  RDF,  cation  - cation  RDF  and  anion  - anion  RDF.  It  is  commonly 
assumed  10  that  the  region  of  the  distances  from  0 through  the  first  maximum  up  to  the  first  minimum  of  the  RDF 
corresponds  to  the  first  neighbourhood  of  the  central  atom,  the  region  from  the  first  minimum  to  the  second  minimum 
corresponds  to  the  second  neighbours,  etc.  The  integrals  of  the  RDF  within  the  limits  given  above  correspond  to  the  number 
of  the  nearest  neighbours  of  the  central  atom,  the  second  nearest  neighbours,  etc. 

Stochastic  Voronoi  polyhedra  seem  to  be  a very  convenient  tool  for  the  description  of  the  long  range  order  in 
disordered  systems  “.  The  definition  of  the  Voronoi  polyhedron  is  illustrated  in  Figure  1.  For  a selected  atom,  referred  to  as 
the  "central  atom",  we  define  the  Voronoi  polygon  as  the  set  of  the  points  that  are  closer  to  the  "central  atom"  than  to  any  of 
the  "medium  atoms"  of  the  disordered  system.  To  find  the  Voronoi  polygon  we  join  the  "central  atom"  and  all  the  "medium 
atoms",  then  at  the  points  halving  the  distances  we  draw  perpendicular  planes  and  select  the  minimal  convex  polyhedron 
constructed  from  these  planes.  Thus,  the  Voronoi  polyhedron  can  be  considered  as  a generalization  of  the  well-known 
Wigner  - Seitz  symmetric  cells  to  the  disordered  systems. 


In  an  ideal  crystal  all  the  Voronoi  polyhedra  are  identical  and  take  on  highly  symmetrical  forms.  When  the  ion 
positions  deviate  from  the  lattice  sites  we  observe  rapidly  increasing  variety  of  the  shapes  of  the  Voronoi  polyhedra.  Some 
of  the  polyhedra  acquire  or  loose  one  or  more  faces.  Thus,  we  can  construct  the  distribution  of  the  number  of  feces  of  the 
Voronoi  polyhedra  drawn  for  the  disordered  set  of  points.  The  sets  that  are  close  to  the  crystalline  order  exhibit  the 
distributions  that  are  sharply  peaked  around  the  number  of  feces  predominating  in  the  system.  In  the  limiting  case  of  the 
perfect  crystal,  the  distribution  is  a Dirac  delta  function  placed  at  fee  number  of  faces  which  corresponds  to  fee  number  of 
fee  first  neighbours  in  fee  crystal  lattice. 

Figure  1 shows  fee  Voronoi  polyhedra  plotted  for  a series  of  cubic  lattices  with  fee  increasing  random  deviations  of 
fee  atom  positions  from  fee  lattice  sites.  The  left  upper  part  shows  fee  polyhedron  for  an  ideal  lattice.  The  upper  right  part  of 
fee  figure  shows  a typical  polyhedron  for  fee  cubic  lattice  when  atoms  are  allowed  20  % deviation  from  fee  lattice  sites.  The 
lower  right  part,  shows  a typical  polyhedron  drawn  for  one  of  fee  configurations  obtained  in  fee  course  of  simulation  of  13 
M LiCl  solution.  It  is  worth  noting  feat  fee  upper  right  and  lower  right  polyhedron  are  very  similar  one  to  another.  These 
polyhedrons  can  be  compared  wife  fee  polyhedron  of  fee  lower  left  part  of  fee  Figure  drawn  for  completely  random  set  of 
points.  We  observed  that  even  minor  deviations  from  fee  ideal  lattice,  of  fee  order  of  1%,  produce  polyhedrons  wife  a wide 
range  of  fee  shapes.  The  distributions  for  fee  distorted  crystals  are  however  clearly  different  than  fee  distribution  for  fee 
random  system. 


It  was  demonstrated  feat  analysis  of  non-sphericity  (or  anisotropy)  factor  of  fee  Voronoi  polyhedra  is  suitable  to  fee 
study  of  the  structure  of  disordered  systems.  The  anisotropy  parameter  is  given  by: 

A3 

a = r-  (2) 

36  n V2 


where  A and  V are  fee  surface  and  volume  of  fee  Voronoi  polyhedron.  For  a sphere  a equals  to  1 , for  fee  bcc,  fee  and  simple 
cubic  lattices  it  equals  to  1.33,  1.35  and  1.91,  respectively.  In  limiting  case  of  the  perfect  crystal,  fee  distribution  of  a is  fee 
delta  Dirac  function.  The  sets  of  points  feat  are  close  to  fee  crystalline  order  exhibit  sharply  peaked  distribution.  For  less- 
ordered  systems  fee  distribution  is  much  broader  and  exhibits  a bell-shaped  form.  The  distributions  of  fee  anisotropy  factors 
of  the  Voronoi  polyhedra  were  calculated  for  fee  NaOH  and  MgCl2  solutions  wife  different  concentrations.  The 
distributions  are  shown  in  Figure  6 together  wife  fee  distribution  for  a random  set  of  points. 
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Fig.l.  Voronoi  polyhedra  plotted  for  a series  of  systems  with  increasing  random  deviations  of  the 
atom  positions  from  the  lattice  sites.  The  left  upper  part:  the  polyhedron  (Wigner-Seitz  cell)  for  an 
ideal  cubic  lattice.  The  upper  right  part:  a typical  polyhedron  for  the  cubic  lattice  when  atoms  are 
allowed  20  % deviation  from  the  lattice  sites.  The  lower  right  part:  shows  a typical  polyhedron  drawn 
for  one  of  the  con  figurati  ons  obtained  in  the  course  of  simulation  of  13  M LiCl  solution.  The  lower  left 
part:  a polyhedron  drawn  for  random  set  of  points. 

3.RESULTS  AND  DISCUSSION 


3.1.  Radial  distribution  functions 

The  ion  - ion  pair  distribution  functions  are  plotted  in  the  following  arrangement:  Figure  2 corresponds  to  NaOH 
solutions  with  concentrations  0.5  M,  3.0  M,  6.0  M and  19.0  M;  Figure  3 to  NaCl  solutions  with  concentrations  0.3  M,  0.9 
M,  2.8  M and  5.1  M;  Figure  4 to  LiCl  solutions  with  concentrations  0.5  M,  5.0  M,  1 1.9  M and  14.0  M;  Figure  5 to  MgCl2 
solutions  with  concentrations  0.2  M,  1.1  M,  2.4  M and  4.9  M.  The  highest  concentrations  of  the  series  are  close  to  the 
saturation  point.  Each  figure  is  composed  of  the  upper,  middle  and  lower  part  representing  the  cation  - anion  distributions, 
the  cation  - cation  distributions  and  the  anion  - anion  distributions,  respectively.  The  distributions  are  additionally  labeled 
in  the  figures. 
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Fig.  4.  The  ion  - ion  radial  distribution  functions  for  LiCl  solutions  with  concentrations:  0.5  M,  5.0  M, 
1 1.9  M and  14.0  M.  The  upper  figure  shows  the  Li+  - Cl'  RDF,  the  figure  in  the  middle  shows  the  Li+  - 
Li+  RDF  and  the  lower  figure  shows  the  Cl"  - Cl'  RDF. 
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The  pair  distribution  function  gives  only  a spherically  averaged  information  about  the  system  structure  and  the 
nonspherical  ionic  associations  are  sometimes  masked.  However,  in  the  case  of  the  concentrated  NaOH,  NaCl  and  LiCl 
solutions  the  presence  of  the  ionic  aggregates  can  be  easily  observed.  We  describe  the  correlations  exhibited  in  Figure  2 (the 
NaOH  case)  in  a more  detailed  way. 

The  cation  - cation  RDF  (Fig.  2,  middle  part)  significantly  changes  with  the  increase  of  the  NaOH  concentration.  In  the 
solution  with  the  lowest  concentration  the  positions  of  sodium  ions  are  not  correlated.  In  the  case  of  3 M solution  there  are 
three  peaks  visible  at  3.2,  6.0  and  9.0  A.  The  first  peak  corresopnds  to  two  Na+  ions  separated  by  a H20  molecule.  The 
position  of  the  first  maximum  moves  to  3.1,  3.0,  2.95  and  2.9  A in  6 M,  10  M,  14  M and  19  M solutions,  respectively. 
Beginning  from  the  6M  NaOH  solution,  the  cation  - cation  RDF’s  exhibit  peak  typical  to  the  crystal-like  arrangement  at  4.2 
A,  corresponding  to  two  sodium  ions  separated  by  two  hydroxyl  ions.  Similar  tendency  could  be  observed  for  the  anion  - 
anion  RDF  plots.  There  are  no  correlations  in  0.5  M solution,  in  3 M (and  more  concentrated)  solutions  the  maxima  at  3.1, 
6.3  and  8.9  A are  appearing  as  the  number  separated  and  Na+  separated  ion  pairs  increases.  The  cation  - anion  pair 
distribution  functions  show  three  peaks  at  2.1  A,  about  5 A and  7.6  A.  The  position  of  the  second  peak  changes  from  5.2  to 
4.7  Awith  the  increase  of  the  concentration.  These  data  suggest  that  with  the  increase  of  the  concentration  the  hydrated  ions 
probably  merge  into  ion-pair  solvated  complexes.  Further  increase  of  the  concentration  leads  to  appearance  of  multi-ionic 
structures  and,  finally,  solvated  crystals.  The  plots  of  the  RDFs  undoubtely  suggest  the  presence  of  higher  ionic  aggregates 
in  concentrated  NaOH  solutions  and,  for  the  highest  concentrations,  the  appearance  of  hydrated  lattice-like  structures. 

Similar  conclusions  result  from  the  analysis  of  the  ion  - ion  radial  distribution  functions  for  NaCl  and  LiCl  solutions 
(Figures  3 and  4,  respectively). 

Figure  5 represents  somewhat  different  case:  magnesium  chloride,  which  is  non-symmetric  solution.  All  the  RDF’s 
change  significantly  with  the  increase  of  the  MgCl2  concentration.  Let  us  begin  with  the  Mg2+  - Mg2,  correlations.  In  the 
solution  with  the  lowest  concentration,  0.2  M,  the  positions  of  magnesium  ions  are  not  correlated.  In  the  case  of  1.1  M 
solution  we  observe  two  peaks  at  relatively  large  distances,  6.8  A and  9.3  A.  From  the  concentration  2.4  M MgCl2  upwards, 
the  cation  - cation  RDFs  exhibit  peaks  typical  to  the  crystal-like  arrangement  at  4.0  A and  4.8  A,  corresponding  to  two 
Mg2+  ions  in  contact  through  the  CL  bridge  and  H20  bridge,  respectively.  Similar  tendency  is  observed  for  the  CL  - CL  RDF. 
We  do  not  observe  any  anion  - anion  correlations  in  0.2  M solution.  From  the  concentration  1.1  M upwards,  the  RDFs 
exhibit  the  maxima  at  3.8  A,  4.9  A which  correspond  to  the  solvent  separated  and  Mg2+  separated  anion  pairs.  The  cation  - 
anion  RDFs  show  three  peaks  at  2.3  A,  around  5.3  A and  7.2  A but  the  peak  at  the  lowest  distance  does  not  appear  in  the  1 . 1 
M solution.  The  structure  of  the  Mg2+  - CL  RDF  at  the  highest  concentrations  resembles  typical  RDFs  for  crystalline  media. 
The  evolution  of  the  RDFs  with  the  increase  of  the  MgCl2  concentration  clearly  reflects  the  process  of  merging  of  the 
solvated  ions  into  solvated  ion-pair  complexes  and  further  into  multi-ion  structures. 

3.2.  Voronoi  polyhedra 

Studies  of  the  radial  distribution  functions  have  been  completed  by  the  calculations  of  the  Voronoi  polyhedra  for  the 
cation  - anion  coordination  in  the  NaOH  and  MgCl2  solutions.  For  a given  cation  (Na+  or  Mg2+)  we  bulid  the  Voronoi 
polyhedron  taking  into  account  all  the  anions  in  the  system.  The  calculation  is  repeated  for  all  the  cations  in  the  simulation 
box  and  then  repeated  for  a large  number  of  the  ion  configurations  in  the  simulation  box  recorded  in  the  course  of  the 
simulation.  The  system  of  the  polyhedra  was  then  analysed  with  the  aim  to  estimate  a degree  of  order  in  the  ionic 
structures.  From  many  possible  versions  of  such  analysis  ( e.g.  plotting  the  distributions  of  the  number  of  walls  of  the 
polyhedra,  calculating  the  effective  coordination  numbers,  etc.  ) we  include  in  the  present  paper  the  distributions  of  the 
anisotropy  factor  of  the  polyhedra  and  their  evolution  with  the  increase  of  the  ion  concentration. 

Figure  6 shows  the  a-parameter  distributions  for  the  NaOH  solutions  with  the  concentrations  0.5  M,  10  M and  19  M 
(upper  part  of  the  figure)  and  for  the  MgCl2  solutions  with  the  concentrations  0.2  M,  2.4  M and  4.9  M (lower  part).  The 
distribution  for  the  random  system  of  the  Lennard  - Jones  particles  is  given  also  in  the  figures.  For  very  low  concentrations 
(not  shown  in  the  figures)  the  a distributions  resemble  the  distributions  for  the  random  system.  Then,  with  the  increase  of 
the  concentration  the  simulated  distributions  steadly  evolve  towards  the  distributions  for  the  distorted  fee  or  bcc  lattices. 
Figure  6 provides  thus  an  additional  support  for  the  presence  of  multi-ion  ordered  structures  in  concentrated  solutions. 
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4.  CONCLUSIONS 


The  results  of  the  simulations  of  the  concentrated  ionic  solutions  allow  us  to  formulate  the  following  conclusion.  For  a 
wide  range  of  of  concentrations,  from  1 M to  20  M,  the  geometry  of  the  ionic  configurations  is  significantly  different  from 
the  geometry  of  a random  set  of  points. 

This  conclusion  is  based  on  the  analysis  of  the  ion  - ion  radial  distribution  functions  and  on  the  comparison  of  the 
characteristics  (such  as  distribution  of  the  anisotropy  coefficent  or  the  distribution  of  the  number  of  faces)  of  the  Voronoi 
polyhedra  constructed  around  the  ions  in  the  simulation  box  to  the  distribution  calculated  for  the  random  tesselation  of 
space.  The  Voronoi  polyhedra  for  the  simulated  ion  configurations  in  solutions  have  similar  characteristics  as  the  polyhedra 
drawn  for  the  slightly  distorted  regular  lattices  of  ions. 

On  a basis  of  these  results,  we  can  put  forward  a conjecture  that  regions  or  domains  of  the  quasi-crystalline  order  of 
ions  exist  in  the  concentrated  solutions.  The  quasi-crystalline  order  is  more  evident  for  lower  part  of  the  concentration 
range,  between  1 M and  5-6  M.  For  the  highest  concentrations,  the  structure  of  the  ionic  configurations  was  determined  not 
only  by  the  ion  - ion  interactions  but,  perhaps  even  predominantly,  by  the  interactions  between  file  ions  and  the  solvent 
molecules  which  belong  at  the  same  time  to  more  than  one  solvation  sphere.  Finally,  the  increase  of  the  charge  of  the  cation, 
like  Mg2*,  results  in  a more  pronounced  quasi-crystalline  order.  It  seems  that  the  results  of  the  simulations  of  the 
concentrated  ionic  solutions  support  to  a large  extent  the  lattice  theory  of  concentrated  electrolytes. 
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